Universal Threshold for the Dynamical Behavior of Lattice Systems with Long-Range 

Interactions 
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Dynamical properties of lattice systems with long-range pair interactions, decaying like l/r" with 
the distance r, are investigated, in particular the time scales governing the relaxation to equilibrium. 
Upon varying the interaction range q, we find evidence for the existence of a threshold at a = d/2, 
dependent on the spatial dimension d, at which the relaxation behavior changes qualitatively and 
the corresponding scaling exponents switch to a different regime. Based on analytical as well as 
numerical observations in systems of vastly differing nature, ranging from quantum to classical, 
from ferromagnetic to antiferromagnetic, and including a variety of lattice structures, we conjecture 
this threshold and some of its characteristic properties to be universal. 
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Most lattice models studied in condensed matter phys- 
ics have interactions of finite range. Despite the long- 
range character of the fundamental electromagnetic in- 
teractions, the presence of positive and negative charges 
gives rise to screening effects that cause interactions to 
be effectively short range and justify an approximation 
by finite-range (and often even nearest-neighbor) inter- 
actions. But a finite-range approximation is not always 
justified. One obvious example are astrophysical systems 
dominated by gravitational interactions where screening 
does not occur. Historically, it was in this context that 
several of the anomalies of long-range interacting systems 
were first discussed, for example negative heat capacities 
of bounded self-gravitating gas spheres [1]. In subsequent 
decades, a general interest in fundamental issues of sta- 
tistical physics of long-range interacting systems arose 
[2]. 

Technically, interactions of finite range turn out to be 
convenient: Exact solutions of many-body models, rare 
as they are, are in most cases restricted to finite-range 
interactions, and also general theorems in thermostatis- 
tical physics (like the convexity properties of thermody- 
namic functions) are often proved under the assumption 
of short-range interactions [3] . "Short-range" here refers 
not only to interactions of finite range, but also to those 
decaying like 1/r" with the distance r and an exponent 
a larger than the spatial dimension d of the system. And 
indeed, various long-range systems, i.e., those with expo- 
nents a < d, have been found to violate ensemble equiv- 
alence, have negative microcanonical specific heat, and 
other peculiarities [2]. 

Out of equilibrium, an intriguing and physically rele- 
vant phenomenon that has been observed in long-range 
systems is the existence of quasistationary states. These 
are nonequilibrium states whose lifetimes t diverge with 
increasing system size N . As a result, for a sufficiently 
large system, relaxation to equilibrium takes place on 



a time scale that is larger than any realistic observa- 
tion time. Such behavior has been reported for various 
long-range systems, ranging from classical toy models 
with mean-field interactions to gravitating systems [4]. 
The exponent a in these studies is usually kept fixed — 
either at a = (corresponding to mean-field interactions) 
where analytical calculations are easier, or at some inte- 
ger number in order to account for classical gravity. 

In this Letter we investigate the relaxation to equi- 
librium of lattice systems with long-range pair interac- 
tions, and in particular the a dependence of the relax- 
ation times. We report analytical as well as numerical 
results on several classes, containing systems of vastly 
differing nature. For all lattice structures studied, and 
regardless of whether the systems are quantum or classi- 
cal, we observe, at a threshold value of a = d/2, a drastic 
change of the relaxational dynamics. In particular, the 
exponent q{a), governing the scaling r ex N'^^"'' of the 
relaxation time r, switches from one regime to another 
at a = d/2. Certain qualitative aspects of the scaling 
laws also appear to be universal, as summarized in Figs. 
1 and 3. Note that the threshold value a = d/2 that is 
relevant for nonequilibrium phenomena differs from the 
one at a = d commonly used for distinguishing between 
long- and short-range behavior in equilibrium statistical 
mechanics. Other aspects of the dynamics are evidently 
nonuniversal, even to the point that relaxation in some 
systems may slow down with increasing system size, but 
accelerate in others. 

This Letter aims at furthering the understanding of 
fundamental aspects of nonequilibrium statistical me- 
chanics of long-range interacting systems. In particu- 
lar, the mechanism of relaxation in long-range interacting 
systems, but also many aspects of its phenomenology, are 
still only poorly understood. Identifying universal prop- 
erties and threshold values may give valuable clues and 
deepen the general understanding. On the more applied 



side, recent developments have made long-range inter- 
acting systems accessible in the laboratory, and an ex- 
perimental verification of some of our results should be 
feasible [5] . Particularly promising for such a check is an 
ion-trap-based quantum simulator as reported by Brit- 
ton et al. [6]. In this setup, a quantum long-range Ising 
model is emulated, and the exponent a governing the 
interaction range can be tuned between and 3. 

Long-range quantum Ising model. — The first class of 
models we are studying consists of spin-1/2 degrees of 
freedom on the N sites of a lattice A, governed by the 
Hamiltonian operator 
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Here, a^ is the z component of the Pauli spin operator, h 
is an external magnetic field, |i— j| denotes the Euclidean 
distance of sites i and j, and the exponent a > deter- 
mines the interaction range. Depending on the sign of the 
coupling J, the spin-spin interactions are ferromagnetic 
(J > 0) or antiferromagnetic ( J < 0). On certain lattices 
(e.g., a triangular lattice), antiferromagnetic interactions 
result in geometrical frustration. An important differ- 
ence to a similar model considered in [7] is the absence of 
an A^-dependcnt normalization factor in front of the first 
sum in (1), introduced in [7] to ensure extcnsivity of the 
energy. We will comment on the role of the factor later 
in this Letter. 

The expectation value (A) of an observable A is given 
by Tr[j4p(t)], where the time evolution of the density 
operator p is governed by the von Neumann equation. 
For all initial density operators po that are diagonal in 
the (T^ tensor product eigenbasis, the time evolution of 
single-spin observables can be computed analytically in 
arbitrary spatial dimension, yielding 
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(see Supplemental Material A). For all finite lattices A, 
(2) is a quasiperiodic function. Proper relaxation to equi- 
librium can be observed in the thermodynamic limit of 
infinite lattice size, where using techniques analogous to 
those in [7-9], we obtain 
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FIG. 1. (Color online) The exponent q of the scaling law r ex 
Af that governs the system-size dependence of the relaxation 
time T. Left: For the long-range quantum Ising model (1) 
the exponent q — min{0, a/d — 1/2} follows from the bounds 
in (3). Right: For the aXY chain (4). The crosses mark data 
points, obtained by a scaling analysis as in Fig. 2. The line is 
plotted as a guide for the eye, indicating two distinct linear 
regimes with a crossover at a/d ~ 0.44. 



the above-mentioned A^-dependent normalization in the 
Hamiltonian (1), which essentially corresponds to taking 
the limits of large A^ and t in a different way. The func- 
tional form of the upper bounds in (3) indeed correctly 
reflects the behavior of |(crf)(i)|, i.e., the powers of A^ 
and t in the exponent agree excellently with a numeri- 
cal evaluation of (2) for large A^, and only the numerical 
constants are, as expected, overestimated (Supplemental 
Material B). 

From Eq. (3) it follows that a change of regime takes 
place at a = d/2. For a < d/2, relaxation to equilibrium 
is Gaussian in time, and the corresponding relaxation 
time T scales like r ex iV"/**"^/^, shrinking to zero in the 
limit of large A^. For a > d/2, relaxation is governed 
by a compressed or stretched exponential in t, with a 
relaxation time that is constant asymptotically for large 
A^. These scaling laws are summarized in Fig. 1 (left). 
The threshold at a = d/2 suggests the following inter- 
pretation: Only for a < d/2 are the pair interactions 
sufficiently long-range such that the relaxation dynamics 
of a single spin is directly influenced, and thereby sped 
up, by the presence of 0{N) spins. Analytical calcula- 
tions of spin-spin correlation functions indicate that fur- 
ther significant qualitative and quantitative changes oc- 
cur in other dynamical quantities [9]. Interestingly, this 
dynamical long-range threshold differs from the equilib- 
rium threshold at a — d below which peculiar long-range 
behavior, like nonequivalence of ensembles or negative 
specific heat, may occur in equilibrium statistical physics. 
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aXY chain. — This model, introduced in [10], consists 
of classical XY spins attached to the sites i of a one- 
dimensional chain and parametrized by the angular vari- 
ables 4>i. The time evolution is generated, via Hamilton's 
equations, by the Hamiltonian function 



valid for large A^ and t (see Supplemental Material B). 
This result differs from the one reported in [7] in a non- 
trivial way. This is a consequence of the absence of 
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where (/> — {4>i, ■ • • , ^a^) is the vector of angle variables 
and p = (pi, . . . ,pn) the vector of conjugate momenta. 
For a = 0, and besides a normalization factor 1/A'' in 
front of the second sum, Eq. (4) reduces to the much stud- 
ied Hamiltonian Mean-Field model [4] , a model known to 
display many of the peculiarities of long-range interact- 
ing systems. In particular, relaxation to equilibrium has 
been studied extensively, and the occurrence of quasista- 
tionary states was observed for large classes of initial con- 
ditions [11]. In equilibrium and for exponents < a < 1, 
the model (4) shows a transition from a magnetized phase 
at energy densities e = E/N < J/4 to an unmagnetized 
phase for e > J/4 [12]. 

Initially we prepare the system in so-called waterbag 
initial distributions, with initial angles (j)i drawn from a 
flat distribution, and initial momenta taking random val- 
ues in the symmetric interval [—A, -|-A] with some A > 0; 
the average energy per particle is then e = A^/6. The 
time evolution is investigated by numerically integrat- 
ing the Hamiltonian equations of motion, using a sixth- 
order symplectic integrator [13]. In earlier studies of the 
a = case, the magnetization m had been monitored 
over time in order to observe relaxation to equilibrium 
[14]. This approach is not viable above the critical en- 
ergy Cc = 1/4 where the initial and the equilibrium value 
of m are both close to zero. For that reason, we moni- 
tor the time evolution of the kurtosis of the momentum 
distribution, k — (p*)/(p^)^, where the angular brackets 
denote averages over the lattice. The kurtosis of the wa- 
terbag initial states we are using is kq — 9/5, whereas 
the Boltzmann equilibrium distribution has Kcq — 3. 

Choosing A — 5/4, we prepare the system at an energy 
density e = 25/96 sa 0.26 slightly above the transition en- 
ergy Cc in the unmagnetized regime. The time evolution 
of the kurtosis k is shown in Fig. 2 (left) for a = 1/8. The 
data reveal a relaxation towards the equilibrium value 
'^cq = 3, on a time scale that depends strongly on the 
system size N. Plotting the same data vs rescaled time 
t/N'^ with q = 1.453, the curves for different N collapse 
onto each other, demonstrating the validity of the scaling 
law; see Fig. 2 (right). Performing such a scaling analysis 
for different values of a, we obtain the scaling exponent 
g as a function of a as shown in Fig. 1 (right). The plot 
reveals two regimes: For < a < ath with ath ~ 0.45, 
q evolves linearly in a as q sa 1.38 -I- 0.5a. The second 
regime, again linear in a, is described by g sa 2.5 — 2a 
for ath < a < 1. This confirms, similar to our findings 
for the quantum Ising model, the presence of two dis- 
tinct power-law regimes for the relaxation times of the 
aXY chain. For a > 1, rescaling of time does not any 
longer lead to a data collapse, so either no such scaling 
law exists in this regime, or larger system sizes would be 
necessary to reach the scaling regime. For other values of 
the energy density e, the data collapse is of lesser quality, 
but the behavior of q{a) is similar to the one shown in 
Fig. 1, though with larger fluctuations (see Supplemen- 
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FIG. 2. (Color online) The kurtosis k, of the momentum distri- 
bution of the aXY chain, plotted for a = 1/8, energy density 
e ^ 0.26, and system sizes N = 256, 512, 1024, 2048 and 4096. 
Data are averaged over 128, 64, 32, 16 and 8 realizations, re- 
spectively. Left: As a function of time t, the relaxation time 
increases with system size A'^. Right: As a function of rescaled 
time t/N'' with scaling exponent q = 1.453. 



tal Material C for simulation data and a more detailed 
discussion) . 

Normalization of the energy scale. — In many papers 
on long-range interacting systems, the pair-interaction 
term in the Hamiltonian is made extensive by equipping 
it with an A'^-dependent normalization factor Af. For the 
models studied in this Letter, such normalization fac- 
tors, as discussed for example in [7] and [10, 12], behave 
as Af ^ N"^'^~^ asymptotically for large N when a < d 
(while no normalization is required for a > d). Exten- 
sivity of the energy is a prerequisite for a well-defined 
thermodynamic limit when making the transition from 
equilibrium statistical mechanics to thermodynamics [3] . 
Moreover, interesting physical phenomena, like equilib- 
rium phase transitions, are caused by competition be- 
tween energetic and entropic effects, and extensivity of 
both these quantities is usually required in order for such 
a competition to survive in the thermodynamic limit. 

For studying the approach to equilibrium, a normaliza- 
tion factor is not necessary, but can be included [15]. An 
A'^-dependent prefactor in the energy induces a further A^ 
dependence of the time scale, shifting the scaling expo- 
nent q of the relaxation time scale to q ~ q + {l — a/d) in 
quantum dynamics, and to g = q-\-(l — a/d)/2 in classical 
Hamiltonian dynamics. The dependence of the exponent 
q on a is shown in Fig. 3 for the long-range quantum 
Ising model and the classical aXY chain. The similari- 
ties between the two models become even more evident 
in this plot, with a crossover from a constant regime for 
a < ath to a linear regime above the threshold. 

In the case of the long-range quantum Ising model, in- 
troducing the normalization factor A^ has the effect of 
even flipping the sign of the exponent, from a negative q 
to a positive q. This implies that, while the relaxation 
to equilibrium takes longer and longer with increasing 
system size in the presence of Af , the opposite happens 
for the original Hamiltonian: Relaxation speeds up with 
increasing A^ , leading to "instantaneous" equilibration in 
the thermodynamic limit. This is different from what is 
observed for the aXY chain, where long-lived quasista- 




and it roughly amounts to requiring that 
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FIG. 3. (Color online) As in Fig. 1, but for scaling exponents 
q as modified by the presence of an A'^-dependent prefactor 
A/" in the Hamiltonian. The left plot is for the long-range 
quantum Ising model, the right for the classical aXY chain. 



tionary states occur independently of whether the factor 
Af is present or not. 

Discussion of the results. — As illustrated in Figs. 1 and 
3, we find that the exponent q in the scaling law r ex N"^ 
of the relaxation time varies linearly as a function of 
a, with a change from one linear regime to another at 
ct = ath- For the quantum Ising model, ath — d/2 is the 
exact location of this change, whereas for the aXY chain 
we find an approximate value of ath ~ 0.45. However, 
results on the largest Lyapunov exponent of the aXY 
model on d-dimensional lattices [16] support the conjec- 
ture that the exact value of ath is d/2 also in the classi- 
cal case: Lyapunov exponents are characteristic quanti- 
ties for the time evolution, quantifying in some sense the 
chaoticity of the dynamics. The authors of [16] find that 
the largest Lyapunov exponent vanishes like N^^, where 
K changes from a constant regime for < a/d < 1/2 to a 
linearly decreasing regime for 1/2 < a/d < 1. The qual- 
itative similarity of K{a) in Fig. 3 of [16] to the plot of 
q{a) in Fig. 3 of the present Letter is striking. At least 
sufficiently far from equilibrium, the sum of positive Lya- 
punov exponents is known to correspond to an entropy 
rate, which in turn provides a link to the speed at which 
equilibrium is approached [17]. When the largest term 
of that sum switches to a different scaling regime, this 
change is expected to reflect also in the sum and, there- 
fore, in the speed at which the system relaxes to equilib- 
rium. These arguments suggest that the transition from 
one linear regime of q{a) to another takes place precisely 
at ath — d/2, not only for the quantum Ising model, but 
also for the aXY chain, and our numerical results are in 
good agreement with this analytical prediction. 

Our understanding of the origin of the universality 
of the threshold at ath — d/2 is partial at best, but 
some physical intuition can be gained from studying Lieb- 
Robinson-type bounds on the propagation of perturba- 
tions. Hastings and Koma [18] report such a bound for a 
broad class of quantum lattice systems with long-range 
interactions. One restriction on the interactions (Eq. 
(2.3) of [18]) is essential for obtaining a nontrivial bound. 



for any given lattice sites i,j G A. By integral approxi- 
mation one finds that (5) is satisfied for a > d/2, repro- 
ducing the threshold value we found for the relaxation 
dynamics. This suggests an appealing, though specula- 
tive, intuitive explanation: In the regime a > ath, re- 
strictions on the speed at which perturbations propagate 
(as given by the Lieb-Robinson bound) are responsible 
for one type of relaxation behavior, whereas the absence 
of such restrictions gives rise to another type of relax- 
ation in the regime a < ath- However, to turn this in- 
tuition into a proof of the universality of the threshold 
at ath, important pieces are still missing, one of them 
being a classical version of the Lieb-Robinson bound for 
long-range interacting systems. 

In summary, substantial analytical as well as numeri- 
cal evidence has been found for the existence of a thresh- 
old at ath — d/2 at which dynamical properties of long- 
range interacting systems show an abrupt change from 
one regime to another. This threshold is found for classi- 
cal Hamiltonian as well as quantum dynamics, on lattices 
of arbitrary dimension and for various lattice structures, 
and for ferromagnetic as well as anti-ferromagnetic in- 
teractions. Anecdotal evidence that the same threshold 
value is of significance also in random models is reported 
in Ref. [19], although here different kinds of randomness 
may behave differently and the scenario becomes more 
involved [20] . Furthermore, a change of regime can be ob- 
served not only in the relaxation times we have studied in 
the present letter, but also in other dynamical quantities, 
be it the above-mentioned largest Lyapunov exponents 
[16] or the emergence of widely separated time scales of 
decaying correlations [9] . Beyond the universality of the 
threshold value ath = d/2, it is tempting to conjecture 
on the basis of Fig. 3 that also the two linear regimes in 
q{a) or q{a) may be model-independent and universal. 
Admittedly, this is speculation, and understanding the 
origin of the observed universality still poses challenges 
for future research. An experimental verification for more 
general quantum Ising models (i.e., in the presence of a 
transverse magnetic field) should be possible with the ion 
trap technology of Ref. [6] . 
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B. Derivation of equation (3) 



A. Derivation of equation (2) 

The derivation is similar to, but simpler than, the one 
reported in Appendix A of [9]. Starting point of the 
calculation is the general expression for the expectation 
value 

(a±)W = Tr(e'^AV±e-'^-Vo) (6) 

with respect to the initial density operator po, where 
we assume po to be diagonal in the ax tensor prod- 
uct eigenbasis. Spin ladder operators are defined as 
CT* = {<T^ ± io'^)/2. Since all terms in the Hamiltonian 
H\ commute, we can factorize the time evolution opera- 
tor, 

exp i-iH^t) = Yl exp (}Jk,i<yl^it) \{ exp {iha^^t) , (7) 

k<l m 

and similarly for the Hermitian conjugate. Here we have 
introduced the notation Jjj- = J /\i — j\°' for the distance- 
dependent couplings. All factors in (7) that do not con- 
tain erf commute with a^ . To compute the time evolu- 
tion in (6), we therefore have to deal with the expression 



Jl exp {-~Uk,^(Tl<7ti) exp {-iha'^t) a^^ 



X exp {iha^t) Y[ exp (iJLiCrfa^t) . (8) 

Making use of [a^'ja^] ~ ±2a^ , the time evolution due 
to the magnetic field h simplifies to 

exp (-iha^t) af exp {iha^t) = a^ exp {^2iht) . (9) 

Picking one lattice site k ^ i, the time evolution of af 
due to the interaction with the spin at k can be written 
as 



exp {-iJk,ialalt) af exp (iJk^talalt) 

= af cos {2tJi^k) + io-f cr| sin (2t J^^fe) 



(10) 



Since the initial state po is assumed to be diagonal in the 
ax tensor product eigenbasis, only diagonal elements of 
the operator e^^'^*a~e~^^^^ in the same basis contribute 
to the trace in (6). For this reason we can drop the second 
term on the right-hand side of (10), as it is proportional 
to a^. Inserting (9) and (10) into (6), we obtain 

{at){t) = (r7±)(0)exp(T2i/ii) J|cos(2tJ,,fe), (11) 

k^i 

where (cf )(0) = Tr (afpo). From (11) we obtain 

«>W = K)(i) + K)(i) 

= (af) (0) cos{2ht) Yl cos (2t J,,fc) . (12) 



This derivation shares some similarities with the ones 
reported in Sec. 4.1 of [7] and Appendix B of [9]. The 
aim is to construct, in the thermodynamic limit of infinite 
lattice size, a nontrivial upper bound on the product 
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occurring in the expression for the modulus of the expec- 
tation value {af){t) in (12). Without loss of generality 
we set I = in the following. 
For any given t, the inequality 
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holds, where 
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and 
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is a d-dimensional ball of radius R, centered around the 
origin. Hence we can split the product (13) in the follow- 
ing way. 
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For finite lattices and some large i, it will happen that 
the remaining product in (17) consists of no factors at 
all, resulting in a trivial upper bound Pa < 1- In fact, 
this was to be expected, since the unitary, quasiperiodic 
dynamics of finite systems gives rise to recurrences. But 
we are interested in the thermodynamic limit where, for a 
given t and large enough lattice A, the remaining product 
in (17) will consist of a large (infinite) number of factors. 
By virtue of the definition of R{t) in (15), the arguments 
of the cosines in these factors are all between —tt/2 and 
7r/2, and for this range of arguments the cosine can be 
bounded by 



|cos(a;)| < 1 -■ 



(18) 



From elementary properties of the exponential function 



it then follows that, in the limit of large lattice size, 
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fceA\BR(t) 

Asymptotically for large lattices of "radius" TZ — N^/'^/2, 
the sum in (19) can be determined by integral approxi- 
mation 
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2^d/2 22a-d^l-2a/d _ (|4Jt/7r|l/" + l^d-2a 

r{d/2) d-2a 

2^(i/2 22a — <i^l — 2Q/d _ 2<i— 2q M j^ /^Id/a — 2 

r((i/2) d-2a 

2^<i/2 J22a-dj^i-2a/d fora<rf/2, 

'^ |d-2a|r((i/2) J2''-2a|4j^/^|d/a-2 fora>rf/2, 

(20) 

where in the last line only the leading contribution in the 
limit of large lattice size N was retained. For the second 
inequality in (20) we have used that 

r^/a ^^\'^^'^" od-2a^d/a-2 

> V, . for x>l. (21) 
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This bound can easily be tightened, at the cost of pushing 
its validity out to larger values of x. Inserting (20) into 
(19) yields the bound 



Pa < 




25 + 2o-d^d/2-2j2 ^^ l_2a/d 

(d-2a)T(d/2) '- ^^ 

2l + d-2o^d/2 |4Jt|d/c 

' {2a-d)T(d/2) I tt I 



for a < d/2, 

for a > d/2, 

(22) 
valid for 4| Ji| > tt and in the limit of large lattice size 
A''. Inserting this bound into (12), we obtain the bound 
on the modulus of the expectation value (erf) as given 
in Eq. (3) of the main paper. A comparison of these 
bounds with an exact evaluation of (13) for finite lattices 
is shown, over more than hundred orders of magnitude, 
in Fig. 4. 

C. More simulation results for the aXY chain 

In the main paper, relaxation time scales of the aXY 
model have been analyzed for e = 25/96 w 0.26, an en- 
ergy density in the unmagnetized phase, located slightly 
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FIG. 4. Comparison of the exact result (13) to the up- 
per bound (22). Left: For a = 1/4 the rescaled logarithm 
Qa = iV^"/'^"MnPA of Pa is plotted as a function of i^ 
(by straightforward evaluation of (13) with Mathematica) . 
With this rescaling, the logarithm of the upper bound in (22) 
is given by the red (upper) straight Hne in the left plot. Qa, 
evaluated for chains of lengths TV = 25001, 50001, 100001, 
and 200001 with open boundary conditions, is shown in blue 
(from top to bottom). All these curves show a linear decay, 
which confirms that the Gaussian decay oc exp(— ct'^) pre- 
dicted by (22) correctly captures the long-time asymptotics 
of Pa, but overestimates, as expected, the constant c. The 
fact that the exact results for Qa do not collapse onto a sin- 
gle curve indicates, however, that the power of A^ in the bound 
(22) is not the sharpest possible one. Right: For a — 3/4, the 
logarithm Qa ~ In Pa (without any A-scaling) is plotted vs. 
rescaled time i"*' ° . With this rescaling, the logarithm of the 
bound in (22) is a straight line (upper red line in the right 
plot). Qa, evaluated for chains of lengths TV = 25001, 50001, 
and 100001 are shown in blue, but the curves fall onto each 
other and are indistinguishable on the scale of the plot. The 
linearly decaying trend in the plot, superimposed by fluctu- 
ations, confirms over a range of more than hundred orders 
of magnitude that the bound in (22) correctly captures the 
compressed or stretched exponential decay of the long-time 
asymptotics of Pa. Again, the numerical constants in the 
compressed or stretched exponential decay are, as expected, 
overestimated. 



above the phase transition energy Ec = 1/4. In the 
present section we will report simulation results also for 
other values of e, both in the magnetized and in the un- 
magnetized phase. The simulation methods used are the 
same as described in the main paper. 



C.l Magnetized phase 

For energy densities e w 0.007 well below e^ — 1/4, 
and starting from unmagnetized out-of-equilibrium ini- 
tial conditions, we observe in our simulations that the 
system experiences a fast transient to a magnetized state, 
and then remains magnetized at all times. This transient 
corresponds to a regime where the macroscopic evolution 
of the system is well captured by a Vlasov equation. The 
initial growth of the magnetization corresponds to a lin- 
ear instability of the Vlasov equation [11], and the expo- 
nential rate neither depends significantly on the system 
size, nor on the interaction range a [21] (see Fig. 5). Con- 
sequently the time needed by the magnetization to reach 
a finite value is determined mainly by the initial condi- 
tion, in particular the initial value of the magnetization 
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FIG. 5. (Color online) The magnetization of the aXY chain, 
plotted for a — 0.0 (thick lines) and a = 0.8 (thin lines), for 
A = 1/5, and system sizes A^ = 256, A^ = 1024, 4096 and 
16384. ty refers to the time of the Vlasov approximation, 
i.e., of the normalized Hamiltonian with the factor A/". 




FIG. 6. (Color online) The kurtosis «£ of the e-distribution 
of the aXY chain, plotted for a = 0.2 and A = 1/5, i.e. 
energy density e ~ 0.007, and system sizes A'^ = 1024, 2048, 
4096, 8192 and 16384. Data are averaged over 5 realizations 
of waterbag initial conditions. Left: As a function of time 
t, the time scale on which Ke relaxes to its equilibrium value 
grows with the system size N. Right: As before, but as a 
function of rescaled time t/N'' with scaling exponent q — 0.6. 



which typically scales as l/y'N. We point out that this 
cannot account for the threshold discussed in the main 
paper. Firstly because, in the mean-field case, the power 
laws of the relaxation times were shown to originate from 
kinetic terms that are not present in the Vlasov equa- 
tion [22]. Secondly, the Vlasov equation holds for the 
Hamiltonian normalized with the factor J\f, so there is 
actually a trivial A^-dependence q — {a — l)/2 for the 
exponential rate, plus the 1/2 from the I/VN of the 
fluctuations; nevertheless, these factors do not yield any 
threshold at a = d/2. 

Apart from the fast transient, the magnetization, tem- 
perature, and kurtosis of the momentum distribution all 
change only very slightly during the relaxation process. 
As a result, their relaxation is virtually impossible to ob- 
serve in the presence of finite-size fluctuations, and we 
need to look for a more suitable observable. To this end 
we introduce the normalized single particle energy at the 
lattice site i, 



J 



N I 

v^ COS 



n_ 



-jI" 



(23) 



whose kurtosis k,^ exhibits a slow relaxation over a large 
range of values (see Fig. 6), which is good basis for ob- 
serving the relaxation dynamics in the presence of fluc- 
tuations. Unfortunately, as can be seen in that same 
figure, the fiuctuations in k^ are also fairly large, and 
this prevents a good data collapse when rescaling time 
as tN~'^^°'\ This suggests that much larger systems need 
to be studied in order to extract scaling laws for the equi- 
libration times in this regime to good accuracy, and this 
will require significantly more computational power. 

Next we study the relaxation time scales at an energy 
density e ~ 0.18 still in the magnetized phase, but only 
slightly below the critical energy. At this energy we ob- 
serve that the system remains unmagnetized for a long 
time before finally switching to the nonzero equilibrium 
value of the magnetization in the magnetized phase (see 
Fig. 7). In this regime, the A^-dependent time scales of 
the Hamiltonian Mean-Field model {aXY model with 
a = 0) were originally studied in [14] and the exponent 
q = 1.2 was extracted by means of a four-parameter fit 
to the magnetization. The result of such a fit depends 
considerably on the choice of the fitting function. 

To avoid the arbitrariness of a fitting function, we use 
a simple rescaling of time as for the other energies stud- 
ied, although it leads to a weaker collapse (see Fig. 7). 
Still, with some optimism, the exponent q of the scaling 
law T ex A^^ can be extracted from the data (see Fig. 
8). Despite the large fluctuations in q, the change of 
regime at around a — d/2 is at least vaguely percepti- 
ble and not inconsistent with the universal threshold we 
propose. The exponent q = 1.2 obtained in [14] for the 
magnetized regime is smaller than the value q — 1.41 
we obtain for a = in the unmagnetized regime, but 
both flndings are consistent with the results of Bouchet 
and Dauxois [22] where — translated into our notation — 
it is proved analytically that q is larger than 1/2 in this 
regime. 




FIG. 7. (Color online) The kurtosis re of the momentum 
distribution of the aXY chain, plotted for a = 0.2 and A = 
1.03, i.e. energy density e ~ 0.18, and system sizes N — 128, 
A^ = 256, N = 512, N = 1024 and 2048. Data are averaged 
over 10 realizations of waterbag initial conditions. Left: As a 
function of time t, the time scale on which k^ relaxes to its 
equilibrium value grows with the system size A''. Right: As 
before, but as a function of rescaled time t/N'' with scaling 
exponent q = 1.3. 
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FIG. 8. (Color online) The exponent q of the scaling law 
r oc A*"* that governs the system-size dependence of the relax- 
ation time T in the aXY long-range chain. The red crosses 
correspond to the unmagnetized regime e ~ 0.26, the blue cir- 
cle to the higher energy regime e « 0.33, the black squares to 
the magnetized regime, for e « 0.007. Data are averaged over 
10 realizations for the e « 0.33 case, and over 5 realizations 
in the e « 0.007 case. 



C.2 Unmagnetized phase 

In the unmagnetized phase above the critical energy, 
we choose the kurtosis n of the momentum distribution 
as an observable, as described in the main paper. In ad- 
dition to the energy value e = 25/96 slightly above the 
phase transition energy Cc = 1/4, we performed simula- 
tions also at higher energy densities, and we observe the 
following antagonistic effects: The higher the energy, the 
smaller are the fluctuations of k in the simulations, and 
the better is the data collapse when rescaling time. How- 
ever, at these higher energies longer simulation times are 
necessary to reach equilibrium. These longer relaxation 
times could be avoided by considering smaller system 
sizes, but this in turn would increase the fluctuations. 
As a consequence, a balance must be found between sim- 
ulation times that can realistically be reached on the one 
hand, and accessible system sizes and/or energies on the 
other hand. The choice e ~ 0.26 as reported in the main 
paper seemed a good compromise and yields a very good 
data collapse of k and small errorbars for the scaling ex- 
ponent q. In Fig. 8 we show results for q{a) at a higher 
energy density e ~ 0.33 (A — 1.4). The results are simi- 
lar to the e « 0.26 case, but with larger fluctuations. 

In summary, for all the energy densities we have sam- 
pled, the simulation results are qualitatively similar to 
the e w 0.26 case discussed in the main paper. The 
results are consistent with the presence of a universal 
threshold value at ath = 1/2 as suggested in the main 
paper, but large fluctuations prevent us from drawing 
stronger conclusions. 
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